rm(list = ls())
load("dataBJPOLS.RData")

inst.1 <- "judgeiv_hd"
endo.1 <- "pti"
outc.1 <- "vote_post"

time.controls <- "as.factor(court_time1) + as.factor(court_time2) + as.factor(court_dow) + as.factor(severity)"
demo.controls <- "age + I(age^2) +  as.factor(race4) + female + vote_pre + as.factor(noteli) + regis_before"
case.controls <- "as.factor(any_drug) +  as.factor(any_weapon) +  as.factor(any_prop) + as.factor(any_prior_case)"

form.3 <- formula(paste(endo.1, "~", time.controls))
form.2 <- formula(paste(inst.1, "~", time.controls))

mod1 <- lm(form.3, data = last.cases)
last.cases$pti_res <- last.cases$pti - predict(mod1, last.cases)

mod1b <- lm(form.2, data = last.cases)
last.cases$judgeiv_res <- last.cases$judgeiv - predict(mod1b, last.cases)

mod2 <- lm(pti_res ~ judgeiv_res, data = last.cases)
judgeiv_res <- seq(min(last.cases$judgeiv_res), max(last.cases$judgeiv_res), by = 0.01)
data.pred2 <- data.table(judgeiv_res)
myPredict <- predict(mod2, newdata = data.pred2, interval = "confidence")

hist(last.cases$judgeiv_res, 
     freq = F, breaks = 10, las = 1,
     xlab = "Leave-out-case Residualized Pretrial Detention Rate",
     main = "", ylim = c(0, 6))
x <- judgeiv_res
ix <- sort(x,index.return=T)$ix

par(new = T)
plot(x, myPredict[, 1], axes=F, xlab=NA, ylab=NA, type = "l", col = "red")
polygon(c(rev(x[ix]), x[ix]), c(rev(myPredict[ ix,3]), myPredict[ ix,2]), col = rgb(0.7,0.7,0.7,0.4) , border = NA)
lines(x[ix], myPredict[ix , 1], col=2, lwd=2)
axis(4, at = seq(from = min(myPredict), to = max(myPredict), by = 0.01), 
     labels = round(seq(from = min(myPredict), to = max(myPredict), by = 0.01), 2), 
     col = 'red', col.axis = 'red', las = 1)
mtext("Residualized Pretrial Detention Rate", side = 4, col = "red", line = 4)

pdf("./Figures/Figure_A04_BJPOLS.pdf", h=8,w=16)    
par(mfrow=c(1,2), mar = c(4, 6, 4, 6), mgp = c(2, 1, 0))

form.3 <- formula(paste(outc.1, "~", time.controls))
mod1 <- lm(form.3, data = last.cases)
last.cases$t_res <- last.cases$vote_post - predict(mod1, last.cases)

mod2 <- lm(t_res ~ judgeiv_res, data = last.cases)

judgeiv_res <- seq(min(last.cases$judgeiv_res), max(last.cases$judgeiv_res), by = 0.01)
data.pred2 <- data.table(judgeiv_res)
myPredict <- predict(mod2, newdata = data.pred2, interval = "confidence")

hist(last.cases$judgeiv_res, 
     freq = F, breaks = 10, las = 1,
     xlab = "Leave-out-case Residualized Pretrial Detention Rate",
     main = "Pretrial Incarceration: Binary Indicator", ylim = c(0, 6), cex.main = 1.5)
x <- judgeiv_res
ix <- sort(x,index.return=T)$ix

par(new = T)
plot(x, myPredict[, 1], axes=F, xlab = NA, ylab = NA, type = "l", col = "red", ylim = c(-0.06, 0.04))
polygon(c(rev(x[ix]), x[ix]), c(rev(myPredict[ ix,3]), myPredict[ ix,2]), col = rgb(0.7,0.7,0.7,0.4) , border = NA)
lines(x[ix], myPredict[ix , 1], col=2, lwd=2)
axis(4, at = seq(from = min(myPredict), to = (max(myPredict)), by = 0.01), 
     labels = round(seq(from = min(myPredict), to = (max(myPredict)), by = 0.01), 2), 
     col = 'red', col.axis = 'red', las = 1)
mtext("Residualized Turnout", side = 4, col = "red", line = 4)

inst.1 <- "judgeiv_cont"
endo.1 <- "log1p(daysjail_t)"
form.3 <- formula(paste(endo.1, "~", time.controls))
form.2 <- formula(paste(inst.1, "~", time.controls))

mod1 <- lm(form.3, data = last.cases)
last.cases$pti_res <- log1p(last.cases$daysjail_t) - predict(mod1, last.cases)

mod1b <- lm(form.2, data = last.cases)
last.cases$judgeiv_res <- last.cases$judgeiv_cont - predict(mod1b, last.cases)

form.3 <- formula(paste(outc.1, "~", time.controls))
mod1 <- lm(form.3, data = last.cases)
last.cases$t_res <- last.cases$vote_post - predict(mod1, last.cases)

mod2 <- lm(t_res ~ judgeiv_res, data = last.cases)

judgeiv_res <- seq(min(last.cases$judgeiv_res), max(last.cases$judgeiv_res), by = 0.01)
data.pred2 <- data.table(judgeiv_res)
myPredict <- predict(mod2, newdata = data.pred2, interval = "confidence")

hist(last.cases$judgeiv_res, 
     freq = F, breaks = 10, las = 1,
     xlab = "Leave-out-case Residualized Pretrial Detention Rate",
     main = "Pretrial Incarceration: Number of Days", ylim = c(0, 1.1), cex.main = 1.5)
x <- judgeiv_res
ix <- sort(x,index.return=T)$ix

par(new = T)
plot(x, myPredict[, 1], axes=F, xlab = NA, ylab = NA, type = "l", col = "red", ylim = c(-0.07, 0.03))
polygon(c(rev(x[ix]), x[ix]), c(rev(myPredict[ ix,3]), myPredict[ ix,2]), col = rgb(0.7,0.7,0.7,0.4) , border = NA)
lines(x[ix], myPredict[ix , 1], col=2, lwd=2)
axis(4, at = seq(from = min(myPredict), to = (max(myPredict)), by = 0.01), 
     labels = round(seq(from = min(myPredict), to = (max(myPredict)), by = 0.01), 2), 
     col = 'red', col.axis = 'red', las = 1)
mtext("Residualized Turnout", side = 4, col = "red", line = 4)

dev.off()
